# Who Do You Trust? - Daniel Goldstein and Johannes Wiedemann
# R-script for Analysis (Main Body)


#====================================================================================================
# load packages

library(tidyverse)
library(readxl)
library(lfe)
library(interactions)
library(estimatr)
library(openxlsx)
library(texreg)
library(dotwhisker)
library(stargazer)
library(cobalt)
library(stringi)



#=======================================================================
# read in data

# Change the working directory as appropriate
setwd("~/../GoldsteinWiedemann_supplementary")



load(file="Data and Data Preparation/GW_2020_MainDataSet_122220.RData") 
load(file="Data and Data Preparation/GW_2020_MRPDataSet_122220.RData")


#=====================================================================================================
#=====================================================================================================





#====================================================================================================
# Analysis starts here




#====================================================================================================
# TABLES (Main Body)


#====================================================================================================


# Create Regression Table for Baseline Diff-in-Diff Effects of Stay-At-Home Orders

panel_1 <- felm(cmi_weekly~post_treatment3|FIPS.y+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_b <- felm(cmi_weekly~post_treatment3+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_c <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_d <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_e <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_1_f <- felm(cmi_weekly~post_treatment3+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


#htmlreg(list(panel_1,panel_1_b,panel_1_c,panel_1_d,panel_1_e,panel_1_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.7,include.rsquared=FALSE,include.groups=FALSE, file = "Tables PoP Final 3/Table1.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE,custom.note = "")

#================================================================


# Create Regression Tables for Stay-at-Home Order Effects by Partisanship

panel_3 <- felm(cmi_weekly~post_treatment3*gop_per_2016|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_b <- felm(cmi_weekly~post_treatment3*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_c <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_d <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_e <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_3_f <- felm(cmi_weekly~post_treatment3*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


#htmlreg(list(panel_3,panel_3_b,panel_3_c,panel_3_d,panel_3_e,panel_3_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Republican Per.","Order x Rep. Per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE, file = "Tables PoP Final 3/Table2.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)

#================================================================




# Create Stay-at-Home Order Effects by Governor Party and Partisanship

panel_14 <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_b <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_c <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_d <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_e <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|week_name|0|State_abb,data = data_weekly_series_1)
panel_14_f <- felm(cmi_weekly~post_treatment3*Republican_gov_2020*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


#htmlreg(list(panel_14,panel_14_b,panel_14_c,panel_14_d,panel_14_e,panel_14_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Republican Governor","Republican per.","Order x Republican Gov.","Order x Rep. per.","Rep. Gov x Rep. per.","Order x Rep. Gov. x Rep. per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.5,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,include.adjrs=FALSE, file = "Tables PoP Final 3/Table3.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)




#================================================================



# Create Stay-at-Home Order Effects by Social Capital

panel_5 <- felm(cmi_weekly~post_treatment3*sk2017_alt|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_b <- felm(cmi_weekly~post_treatment3*sk2017_alt+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_c <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_d <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_e <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_5_f <- felm(cmi_weekly~post_treatment3*sk2017_alt+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


#htmlreg(list(panel_5,panel_5_b,panel_5_c,panel_5_d,panel_5_e,panel_5_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Social Capital","Order x Soc. Cap.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE, file = "Tables PoP Final 3/Table4.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)






#================================================================



# Create Regression Tables for Stay-at-Home Order Effects by Social Capital and Partisanship

panel_15 <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_b <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_c <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_d <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_e <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_15_f <- felm(cmi_weekly~post_treatment3*sk2017_alt*gop_per_2016+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])


#htmlreg(list(panel_15,panel_15_b,panel_15_c,panel_15_d,panel_15_e,panel_15_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Social Capital","Republican per.","Order x Soc. Cap.","Order x Rep.","Soc. Cap. x Rep.","Order x Soc. Cap. x Rep.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.5,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,include.adjrs=FALSE,file = "Tables PoP Final 3/Table5.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)





#================================================================


# Create Regression Tables for Stay-at-Home Order Effects by Vaccination Rate

panel_20 <- felm(cmi_weekly~post_treatment3*analysis_value|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_b <- felm(cmi_weekly~post_treatment3*analysis_value+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_c <- felm(cmi_weekly~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_d <- felm(cmi_weekly~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_e <- felm(cmi_weekly~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices+cmi_weekly_2019|State_abb+week_name|0|State_abb,data = data_weekly_series_1)
panel_20_f <- felm(cmi_weekly~post_treatment3*analysis_value+Age65AndOlderPct2010+log(PopDensity2010+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct+WhiteNonHispanicPct2010+BlackNonHispanicPct2010+AsianNonHispanicPct2010+RuralUrbanContinuumCode2013+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|State_abb+week_name|0|State_abb,weights = data_weekly_series_1$TotalPopEst2018[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE],data = data_weekly_series_1[is.na(data_weekly_series_1$TotalPopEst2018)==FALSE,])



#htmlreg(list(panel_20,panel_20_b,panel_20_c,panel_20_d,panel_20_e,panel_20_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Order","Vaccination Per.","Order x Vacc. Per.","Income p.c.","Unemployment","Education","Rur.-Urb.","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian","CMI 2019"),caption = "",caption.above = TRUE,scalebox = 0.65,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,file = "Tables PoP Final 3/Table6.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)



#================================================================


# Create Regression Tables for Stay-at-Home Order Effects by Trust in Others (MRP)


panel_0 <- felm(cmi_new_weekly~TrustOthers|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1 <- felm(cmi_new_weekly~post_treatment3*TrustOthers|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_b <- felm(cmi_new_weekly~post_treatment3*TrustOthers+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_c <- felm(cmi_new_weekly~post_treatment3*TrustOthers+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_d <- felm(cmi_new_weekly~post_treatment3*TrustOthers+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_f <- felm(cmi_new_weekly~post_treatment3*TrustOthers+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,weights = MRP_weekly_series_1$TotalPopEst2018.x[is.na(MRP_weekly_series_1$TotalPopEst2018.x)==FALSE],data = MRP_weekly_series_1[is.na(MRP_weekly_series_1$TotalPopEst2018.x)==FALSE,])


#htmlreg(list(panel_0,panel_1,panel_1_b,panel_1_c,panel_1_d,panel_1_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Soc Trust","Order","Order x Soc Trust","Income p.c.","Unemployment","Education","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian"),caption = "",caption.above = TRUE,scalebox = 0.7,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,file = "Tables PoP Final 3/Table7.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)


#================================================================


# Create Regression Tables for Stay-at-Home Order Effects by Political Trust (MRP)

panel_0 <- felm(cmi_new_weekly~NoCorruption_num|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1 <- felm(cmi_new_weekly~post_treatment3*NoCorruption_num|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_b <- felm(cmi_new_weekly~post_treatment3*NoCorruption_num+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_c <- felm(cmi_new_weekly~post_treatment3*NoCorruption_num+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_d <- felm(cmi_new_weekly~post_treatment3*NoCorruption_num+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,data = MRP_weekly_series_1)
panel_1_f <- felm(cmi_new_weekly~post_treatment3*NoCorruption_num+Age65AndOlderPct2010.x+log(PopDensity2010.x+1)+log(PerCapitaInc)+UnempRate2018+Ed5CollegePlusPct.x+WhiteNonHispanicPct2010.x+BlackNonHispanicPct2010.x+AsianNonHispanicPct2010.x+PctEmpAgriculture+PctEmpManufacturing+PctEmpServices|week_name|0|State_abb,weights = MRP_weekly_series_1$TotalPopEst2018.x[is.na(MRP_weekly_series_1$TotalPopEst2018.x)==FALSE],data = MRP_weekly_series_1[is.na(MRP_weekly_series_1$TotalPopEst2018.x)==FALSE,])

#htmlreg(list(panel_0,panel_1,panel_1_b,panel_1_c,panel_1_d,panel_1_f),stars = c(0.01,0.05,0.1),custom.coef.names = c("Pol Trust","Order","Order x Pol Trust","Income p.c.","Unemployment","Education","Pct Agr.","Pct Manuf.","Pct Serv.","Pct over 65","Pop. Dens.","Pct White","Pct Black","Pct Asian"),caption = "",caption.above = TRUE,scalebox = 0.7,include.rsquared=FALSE,custom.note = "",include.groups=FALSE,file = "Tables PoP Final 3/Table8.doc", 
 #       inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
  #      head.tag = TRUE, body.tag = TRUE)





#=============================================================================
#=============================================================================
# FIGURES (Main Body)


#================================================================

# Create heat map that display the mobility change in week 15 of 2020 and 2019
library(usmap)

dummy <- data_weekly_series_1 %>% filter(week_name=="2020-W15") %>% mutate(fips=as.numeric(FIPS.y),cmi_change=(cmi_weekly-cmi_weekly_2019)) %>% select(fips,cmi_change) %>% ungroup() 
plot_usmap(regions = "counties",data = dummy,values="cmi_change", color=NA, size=0.00001)+ 
  scale_fill_continuous(
    low = "red", high = "white", name = "CMI Change", label = scales::comma
  ) + theme(legend.position = "right")


#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure1_HeatMap1.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')



#================================================================


# Create heat map that illustrates when states implemented stay-at-home orders
library(usmap)

dummy <- data_weekly_series_1 %>%distinct(State_abb,.keep_all = TRUE) %>% mutate(state=State_abb) %>% select(state,TreatmentWeek_num)

plot_usmap(data = dummy, values = "TreatmentWeek_num", labels=FALSE) + 
  scale_fill_continuous( low = "white", high = "orange", 
                         name = "Shelter Order Week", label = scales::comma
  ) + 
  theme(legend.position = "right") + 
  theme(panel.background = element_rect(colour = "white")) 

#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure2_HeatMap.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')




#================================================================


# Create Parallel Trends for Mobility Trends for Selected States

data_weekly_series_1 %>%  mutate(cmi_counterf=ifelse(TreatmentWeek_rel2<0,cmi_weekly,NA))%>% ungroup() %>%group_by(week_num)%>% mutate(cmi_counterf2=mean(cmi_counterf,na.rm=TRUE))%>%filter(NeverTreated_CA_MI_NC!="Rest") %>%mutate(cmi_weekly_pre=ifelse(week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>% mutate(cmi_weekly_post=ifelse(week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%mutate(ExampleStates=as.factor(NeverTreated_CA_MI_NC))%>% 
  mutate(cmi_counterf2=ifelse(post_treatment3==0,(sum(cmi_counterf,na.rm=TRUE)-cmi_counterf)/(n()-1),cmi_counterf2))%>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=cmi_weekly_pre,linetype="solid",color=ExampleStates),geom = 'line',size=1.1) +
  stat_summary(aes(y=cmi_weekly_post,linetype="dashed",color=ExampleStates),geom='line',size=1.1) +
  theme_minimal()+
  xlim(8,16)+
  xlab("Week Number")+
  ylab("Mobility")+
  scale_color_manual(name="Example States",breaks = c("CA", "MI", "NC","Never"),
                     values=c("red", "green", "blue","grey"))+
  scale_linetype_discrete(name="Treatment Status",labels=c("Post","Pre"))

#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure3_ParallelTrendsSelectedStates.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')



#================================================================


# Create Parallel Trends Figure for Mobility Trends by Treatment Week

plot1 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_12_pre=ifelse(TreatmentWeek_num2==12&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_12_post=ifelse(TreatmentWeek_num2==12&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(TreatedDummy=if_else(TreatmentWeek_num2==12,1,0))%>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=group_12_pre,linetype="solid"),geom = 'line',color="red") +
  stat_summary(aes(y=group_12_post,linetype="dashed"),geom='line',color="red") +
  stat_summary(aes(y=pre_3,linetype="dotted"),geom='line',color="grey") +
  stat_summary(aes(y=pre_2,linetype="longdash"),geom='line',color="grey") +
  geom_vline(aes(xintercept = 11),color="black",linetype="dashed")+
  theme_minimal()+
  theme(legend.position = "bottom")+
  guides(size=FALSE,color=FALSE)+
  scale_linetype_discrete(name="Status",labels=c("Pre","Post","Late","Never"),breaks=c("solid","dashed","longdash","dotted"))+
  xlim(8,16)+
  xlab("Week")+
  ylab("Mobility")+
  ggtitle("Week 12")

plot2 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_13_pre=ifelse(TreatmentWeek_num2==13&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_13_post=ifelse(TreatmentWeek_num2==13&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=group_13_pre,linetype="solid"),geom = 'line',color="red") +
  stat_summary(aes(y=group_13_post,linetype="dashed"),geom='line',color="red") +
  stat_summary(aes(y=pre_3,linetype="dotted"),geom='line',color="grey") +
  stat_summary(aes(y=pre_2,linetype="longdash"),geom='line',color="grey") +
  geom_vline(aes(xintercept = 12),color="black",linetype="dashed")+
  theme_minimal()+
  guides(size=FALSE,color=FALSE)+
  theme(legend.position = "bottom")+
  scale_linetype_discrete(name="Status",labels=c("Pre","Post","Late","Never"),breaks=c("solid","dashed","longdash","dotted"))+
  xlim(8,16)+
  xlab("Week")+
  ylab("Mobility")+
  ggtitle("Week 13")

plot3 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_14_pre=ifelse(TreatmentWeek_num2==14&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_14_post=ifelse(TreatmentWeek_num2==14&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(TreatmentWeek_num2==15&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  mutate(Treated=ifelse(TreatmentWeek_num2==14,1,0)) %>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=group_14_pre,linetype="solid"),geom = 'line',color="red") +
  stat_summary(aes(y=group_14_post,linetype="dashed"),geom='line',color="red") +
  stat_summary(aes(y=pre_3,linetype="dotted"),geom='line',color="grey") +
  stat_summary(aes(y=pre_2,linetype="longdash"),geom='line',color="grey") +
  geom_vline(aes(xintercept = 13),color="black",linetype="dashed")+
  theme_minimal()+
  guides(size=FALSE,color=FALSE)+
  theme(legend.position = "bottom")+
  scale_linetype_discrete(name="Status",labels=c("Pre","Post","Late","Never"),breaks=c("solid","dashed","longdash","dotted"))+
  xlim(8,16)+
  xlab("Week")+
  ylab("Mobility")+
  ggtitle("Week 14")


plot4 <- data_weekly_series_1 %>% mutate(early_late_never=ifelse(TreatmentWeek_num2%in% c(12,13),1,ifelse(TreatmentWeek_num2%in%c(14,15),2,3)))%>% 
  mutate(group_15_pre=ifelse(TreatmentWeek_num2==15&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),group_15_post=ifelse(TreatmentWeek_num2==15&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(pre_1=ifelse(early_late_never==1&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_2=ifelse(early_late_never==2&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA),pre_3=ifelse(early_late_never==3&week_num<=(TreatmentWeek_num2-1),cmi_weekly,NA))%>%
  mutate(post_1=ifelse(early_late_never==1&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_2=ifelse(early_late_never==2&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA),post_3=ifelse(early_late_never==3&week_num>=(TreatmentWeek_num2-1),cmi_weekly,NA)) %>%
  ggplot(aes(x=week_num)) +
  stat_summary(aes(y=group_15_pre,linetype="solid"),geom = 'line',color="red") +
  stat_summary(aes(y=group_15_post,linetype="dashed"),geom='line',color="red") +
  stat_summary(aes(y=pre_3,linetype="longdash"),geom='line',color="grey") +
  geom_vline(aes(xintercept = 14),color="black",linetype="dashed")+
  theme_minimal()+
  theme(legend.position = "bottom")+
  scale_linetype_discrete(name="Status",labels=c("Pre","Post","Never"),breaks=c("solid","dashed","longdash"))+
  guides(size=FALSE,color=FALSE)+
  xlim(8,16)+
  xlab("Week")+
  ylab("Mobility")+
  ggtitle("Week 15")


library(cowplot)
plot_grid(plot1, plot2,plot3,plot4, ncol=2)





#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure4_ParallelTrendsbyTreatmentWeek.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')


#================================================================


# Produce mobility data by partisanship plots (for each state separately)

dummy <- data_weekly_series_1 %>% filter(is.na(State_abb)==FALSE) %>% filter(week_num>8)%>% filter(is.na(gop_per_2016)==FALSE)%>% mutate(cmi_weekly_Rep=ifelse(gop_win_2016==1,cmi_weekly,NA))%>%mutate(cmi_weekly_Dem=ifelse(gop_win_2016==0,cmi_weekly,NA)) %>%select(State_abb,FIPS.y,cmi_weekly,week_num,gop_per_2016,gop_win_2016,cmi_weekly_Rep,cmi_weekly_Dem,TreatmentWeek_num) %>% ungroup() 
ggplot(data = dummy, aes(week_num, cmi_weekly)) +
  theme_bw()+
  theme(legend.position="bottom",legend.box="vertical",element_blank())+
  geom_jitter(aes(color=gop_per_2016),size=1)+
  labs(
       y = "Mobility", x = "Week") + 
  geom_smooth(aes(y=ifelse(week_num<=(TreatmentWeek_num-1),cmi_weekly_Rep,NA),x=week_num),method = "lm",se=FALSE,color="tomato1")+
  geom_smooth(aes(y=ifelse(week_num>=(TreatmentWeek_num-1),cmi_weekly_Rep,NA),x=week_num),method = "lm",se=FALSE,color="tomato1")+
  geom_smooth(aes(y=ifelse(week_num<=(TreatmentWeek_num-1),cmi_weekly_Dem,NA),x=week_num),method = "lm",se=FALSE,color="navyblue")+
  geom_smooth(aes(y=ifelse(week_num>=(TreatmentWeek_num-1),cmi_weekly_Dem,NA),x=week_num),method = "lm",se=FALSE,color="navyblue")+
  geom_smooth(aes(y=ifelse(is.na(TreatmentWeek_num),cmi_weekly_Rep,NA),x=week_num),method = "lm",se=FALSE,color="tomato1")+
  geom_smooth(aes(y=ifelse(is.na(TreatmentWeek_num),cmi_weekly_Dem,NA),x=week_num),method = "lm",se=FALSE,color="navyblue")+
  geom_vline(aes(xintercept = (TreatmentWeek_num-1)),color="black",linetype="dashed")+
  scale_color_gradient(
    low = "navyblue", high = "tomato1", name = "Rep. Vote Share")+
  facet_wrap(~ State_abb)

#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure5_PartisanPlotsFacets.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')





#================================================================


# Create Mobility Trends by Partisanship Figures

data_weekly_series_1 %>% filter(is.na(gop_win_2016)==FALSE)%>% ggplot(aes(TreatmentWeek_rel2, cmi_weekly,color=as.factor(gop_win_2016))) +
  stat_summary(geom = 'line',size=1.1) +
  geom_vline(xintercept = -1,color="purple",linetype="dashed") +
  theme_minimal()+
  xlim(-5,4)+
  xlab("Week Number Relative to Treatment")+
  ylab("Mobility")+
  scale_color_manual(breaks = c("1", "0"), values=c("navyblue", "tomato1"),labels = c("Republican", "Democratic"),name="Party")

#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure6_MobilityTrendsByPartisanship.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')




#=====================================================================================================

# Produce mobility data by social capital plots (for each state separately)



dummy <- data_weekly_series_1 %>% filter(is.na(State_abb)==FALSE) %>% filter(week_num>8)%>% filter(is.na(sk2017_alt)==FALSE)%>% mutate(cmi_weekly_High=ifelse(sk2017_alt_bin==1,cmi_weekly,NA))%>%mutate(cmi_weekly_Low=ifelse(sk2017_alt_bin==0,cmi_weekly,NA)) %>%select(State_abb,FIPS.y,cmi_weekly,week_num,sk2017_alt,sk2017_alt_bin,cmi_weekly_High,cmi_weekly_Low,TreatmentWeek_num) %>% ungroup() 
ggplot(data = dummy, aes(week_num, cmi_weekly)) +
  theme_bw()+
  theme(legend.position="bottom",legend.box="vertical",element_blank())+
  geom_jitter(aes(color=sk2017_alt),size=1)+
  labs(
       y = "Mobility", x = "Week") + 
  geom_smooth(aes(y=ifelse(week_num<=(TreatmentWeek_num-1),cmi_weekly_High,NA),x=week_num),method = "lm",se=FALSE,color="darkgreen")+
  geom_smooth(aes(y=ifelse(week_num>=(TreatmentWeek_num-1),cmi_weekly_High,NA),x=week_num),method = "lm",se=FALSE,color="darkgreen")+
  geom_smooth(aes(y=ifelse(week_num>=(TreatmentWeek_num-1),cmi_weekly_Low,NA),x=week_num),method = "lm",se=FALSE,color="gold")+
  geom_smooth(aes(y=ifelse(week_num<=(TreatmentWeek_num-1),cmi_weekly_Low,NA),x=week_num),method = "lm",se=FALSE,color="gold")+
  geom_smooth(aes(y=ifelse(is.na(TreatmentWeek_num)==TRUE,cmi_weekly_Low,NA),x=week_num),method = "lm",se=FALSE,color="gold")+
  geom_smooth(aes(y=ifelse(is.na(TreatmentWeek_num)==TRUE,cmi_weekly_High,NA),x=week_num),method = "lm",se=FALSE,color="darkgreen")+
  geom_vline(aes(xintercept = (TreatmentWeek_num-1)),color="black",linetype="dashed")+
  scale_color_gradient(
    low = "gold", high = "darkgreen", name = "Social Capital")+
  facet_wrap(~ State_abb)


#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure7_SocCapPlotsFacets.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')


#================================================================


# Create Mobility Trends by Social Capital Figure

data_weekly_series_1 %>% filter(is.na(sk2017_alt_bin)==FALSE)%>% ggplot(aes(TreatmentWeek_rel2, cmi_weekly,color=as.factor(sk2017_alt_bin))) +
  stat_summary(geom = 'line',size=1.1) +
  geom_vline(xintercept = -1,color="purple",linetype="dashed") +
  theme_minimal()+
  scale_color_manual(breaks = c("1", "0"), values=c("darkgreen", "gold"),labels = c("High", "Low"),name="Social Capital")+
  xlim(-5,4)+
  xlab("Week Number Relative to Treatment")+
  ylab("Mobility")

#ggsave("Figures PoP Final 3/GoldsteinWiedemann_Figure8_MobilityTrendsBySocialCapital.tiff",height = 6,width = 8.25,units = "in",dpi = 1000,device='tiff')



#====================================================================================================
#=====================================================================================================